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Abstract 



^^ We investigate how scale-free (SF) and Erdos-Renyi (ER) topologies affect the interplay between evolvability and robustness of 
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model gene regulatory networks with Boolean threshold dynamics. In agreement with |Oikonomou and Cluzel (2006) we find that 
networks with SFi,, topologies, that is SF topology for incoming nodes and ER topology for outgoing nodes, are significantly more 
evolvable towards specific oscillatory targets than networks with ER topology for both incoming and outgoing nodes. Similar results 
are found for networks with SFboth and SFout topologies. The functionality of the SFout topology, which most closely resembles the 
structure of biological gene networks (iBab u et aL] |2004| l, is compared to the ER topology in further detail through an extension to 
multiple target outputs, with either an oscillatory or a non-oscillatory nature. For multiple oscillatory targets of the same length, 
the differences between SFout and ER networks are enhanced, but for non-oscillatory targets both types of networks show fairly 
similar evolvability. We find that SF networks generate oscillations much more easily than ER networks do, and this may explain 
why SF networks are more evolvable than ER networks are for oscillatory phenotypes. In spite of their greater evolvability, we 
find that networks with SFom topologies are also more robust to mutations than ER networks. Furthermore, the SFout topologies 
are more robust to changes in initial conditions (environmental robustness). For both topologies, we find that once a population of 
networks has reached the target state, further neutral evolution can lead to an increase in both the mutational robustness and the 
environmental robustness to changes in initial conditions. 

Key words: gene regulatory networks. Boolean threshold networks, robustness, evolvability, scale-free, oscillatory gene 
expression 
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1. Introduction 

1.1. Robustness &■ evolvability in genetic regulatory networks 

Biological systems require their phenotype to be robust to 
a variety of perturbations. They must be mutationally robust 
to minimize the possibility of a permanent deleterious muta- 
tion in the genome. They must be environmentally robust to 
mitigate the effects of temporary macroscopic environmental 
variation, or stochastic noise in microscopic biochemical pro- 
cesses. By possessing a robust phenotype, functions required 
for survival can be performed consistently and descendants can 
expect a similarly reliable genotype (Wagner 2005 | l. 

At the same time, environments can change permanently and 
populations must adapt by producing heritable phenotypic vari- 
ation. This ability to generate phenotypic innovation from ge- 
netic changes under natural selection is often called evolvabil- 



ity, and is observed throughout the biological world ( Dawkins 



T989||Kirschner and Gerhart,,1998; Wagner||2005] l. It may it- 
self be an evolvable trait ( Earl and Deem 2004| i. The relation- 
ship between environmental robustness and evolvability is very 
complex ( West-Eberhard 2003; Kirschner and Gerhartl 1200511. 
On the other hand, it would seem at first sight that evolvability 
and mutational robustness are straightforwardly antagonistic: 



the more easily mutations change the properties of an organ- 
ism, the less mutationally robust it should be. 

However, this simple picture can be challenged on numerous 



fronts (|Bloom et al.||2007a|b| |Wagner| |2008b|a| [Daniels et alT 



|2008[[Draghi et al. 2010 1. It is important to take into account 
the fact that evolution acts on populations and not on individ- 
uals. Furthermore, the structure of the evolutionary landscape, 
the population size Npop, the mutation rate fi and the type of 
phenotypic variation required are all factors that can influence 
the relationship between mutational robustness and evolvabil- 
ity. 

For example Wagner, using RNA secondary structure as a 
model (Wagner, 2008b), demonstrated the importance of dis- 
tinguishing between genotype (sequence) robustness and evolv- 
ability, which do share an antagonistic relationship, and pheno- 
type (structure) robustness and evolvability which, by contrast, 
share a beneficial relationship. For a single sequence (geno- 
type), mutational robustness is simply the fraction of mutations 
that alter the phenotype. However, a single RNA phenotype (a 
particular secondary structure) may be the minimal free energy 
structure for a large number of different genotypes. Such a set 
of genotypes is called a neutral space, and during evolution, a 
population can spread on the neutral space to generate a large 
amount of genetic diversity. The larger the neutral space, the 
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more likely that random mutations will generate another mem- 
ber of the neutral space, so the larger the mutational robustness. 



It can be demonstrated (Wagner 2008b) that for RNA, larger 



neutral spaces have a larger number of different phenotypes ac- 
cessible with a single mutation. Thus a population that can 
generate greater genetic diversity through non-deleterious mu- 
tations across a larger neutral space, also has more evolutionary 
innovation available, making such phenotypes more evolvable. 
Furthermore, for RNA, the utility of the neutral space with re- 
spect to robustness and evolvability is shown to be independent 
of the evolutionary dynamical regime (being of importance in 
both the Npopfi » 1 and Npopiu <K 1 cases, where Npop is popu- 
lation size and fi is mutation rate). 

These arguments depend on the structure of the neutral 



space (Wagner 2008a I. If it is made up of many disconnected 
pockets, then a population may not be able to explore the en- 
tire neutral space. If on the other hand all the neutral neigh- 
bours are clustered together tightly in one part of the space, 
then they may not have a large diversity of other phenotypes 
within easy reach. Models of RNA secondary structure have the 
advantage that they can be relatively easily analyzed by com- 



putational methods (Wagner 2005 i. Thus detailed questions 
about the topology of neutral spaces and the effects of changing 
mutation rates and population sizes can be analyzed in some 
detail. For other systems, however, the relationship between 
mutational robustness, neutral spaces and evolvability may be 



more complex (Draghi et al. 2010[ ) 



In this paper we focus on gene regulatory networks (GRNs), 
another set of biological systems that can be analyzed compu- 
tationally. There has been an enormous amount of interest in 
modeling GRNs because they play a central role in biological 
functionality ( |Alon[ |2006| l. Moreover, it is increasingly being 
recognized that altering gene expression patterns is a highly ef- 
fective evolutionary mechanism. Important phenotypic changes 
are often achievable through very small numbers of mutations 
to regulatory regions (Carroll 20051. For these reasons, the 
mutational robustness, environmental robustness and evolvabil- 
ity of GRNs have received much theoretical attention. For 
example, the reliability under perturbed system configurations 
(environmental robustness) has been investigated in (Aldana , 
[20031 iKlemm and Bomhoidt|[2005l l whilst the evolution of this 



type of robustness is considered in Braunewell and Bornholdt 
( |2009) . Robustness to mutation has also been shown to be an 
evolvable property, correlating with robustness to initial gene 
expression states (Ciliberti et al. 2007b). Furthermore, there is 



evidence suggesting that neutral spaces can enhance evolvabil- 
ity of expression patterns (Ciliberti et al. 2007a| l. The effect of 
varying environments has been demonstrated to increase evolv- 



ability ( Crombach and Hogeweg 2008) , whilst the role of over- 
all network topology has also been considered (|Aldana[ |2003 



[Oikonomou and Cluzel| [20061 |Aldana"eraLl [2007^ 



1.2. Modelling GRNs 

Realistically modelling large complex interacting GRNs is a 
difficult task. Even a modest sized network can have an enor- 
mous number of parameters depending on which details of the 
fundamental processes involved in transcription and translation 



are taken into account (Alon 2006| l. In practice, therefore, 
many different approaches to simplify the dynamics of GRNs 
can be found in the literature, each with their strengths and lim- 
itations ( Polynikis et al.[[2009| l. 

A particularly simpUfied model, the Boolean network, was 
proposed by Kauf fman[ ( |1969| l. The generic Boolean model 
makes two important assumptions. First, the gene states are the 
variables under time evolution, as opposed to the concentrations 
of gene products. These states are further assumed to be either 
active or inactive, i.e. either "on" or "off". This approximation 
is based upon evidence that the changes in gene state are co- 
operative transitions, allowing the approximation of sigmoidal 



functions as step- wise ones (Kauffman 1993) . The second as- 
sumption is synchronous updating of gene states throughout the 
network in discrete time steps through physically reasonable 
rules. In the original model, these rules were random Boolean 
functions of all incoming nodes, whose output determined a 
given node's state in the next time step. Other update rules, 
based upon neural threshold network models ([Hopfield 1982|l, 



have also been applied in Boolean threshold networks (Wagner 



[1994| l. The Boolean model has been extensively used in theo- 
retical biology and, despite its coarse-grained nature, has been 
successfully applied, for example, to the yeast cell cycle net- 



work (Li et al.l 2004), the segment polarity gene network of D. 



melanogaster (Albert and Othmer,,2003) and the floral cell fate 
of A. thaliana ( |Espinosa-Soto et aL 2004| l. 

An important feature of Boolean network models is the man- 
ifestation of attractors. Following a certain amount of time 
evolution of the system, all initial configurations converge on 
a very small subset of possible configurations - the attractors 
of the network, also referred to as the network's attractor land- 
scape ( |KauffTnan[ [I993| l. Boolean network models have been 
observed to have ordered, critical and chaotic phases correlat- 
ing with an increase in average connectivity (Kauffman 1993| l. 
In the ordered phase, small perturbations in an initial configu- 
ration tend not to spread, with the same attractor reached as in 
the unperturbed case (Aldana[ [2003| l. However, in the chaotic 
phase these small perturbations spread, possibly resulting in the 
network arriving at a different attractor The critical phase lies 
between the ordered and chaotic phases. 

It has been hypothesised that the attractor landscape rep- 
resents the possible "cell fates" given a cell's regulatory in- 
teractions. Different sets of initial configurations induce the 
production of different attractors, regarded as different cell 
types ([Kauffman 1993|l. Recently, some evidence has been 



produced supporting this theory (Espinosa-Soto et al. 2004 



[Huang et all 2005 | l. Within this interpretation, there is a diverse 
range of possibilities, from "frozen" attractors of unit length, 
as required in developmental gene expression patterns, to ones 
including longer oscillatory periods capable of more complex 
function, for example, the changes in gene expression over the 



course of the cell cycle (Li et al. 2004 1. 



1.3. The role of topology 

Within the Boolean approximation that coarse-grains out 
many biochemical parameters, a GRN can be viewed as a di- 
rected graph, with genes as nodes and gene interactions as 



Table 1 : The topologies investigated are defined with respect to combinations of 
incoming and outgoing node degree distributions. When average connectivity 
is small compared with the number of nodes, the binomial distribution may be 
approximated as a Poisson distribution. 



Topology 


Incoming degree 


Outgoing degree 


ER 


Binomial/Poisson 


Binomial/Poisson 


^i^out 


Binomial/Poisson 


Power-law 


SFi„ 


Power-law 


Binomial/Poisson 


SFboth 


Power-law 


Power-law 



edges. The number of incoming and outgoing connections for 
a given node — its incoming/outgoing degree — can be de- 
scribed statistically by distributions that characterize the overall 



structures or topologies (Barabasi and Oltvai 2004 1. If a set of 



nodes is randomly assigned connections, the resulting degrees 



are binomially distributed, forming Erdos-Renyi graphs ( Erdos 



and Renyi I960i. Another class is the scale-free (SF) topol- 
ogy, where the node degrees are power-law distributed. Given 
that the incoming and outgoing degrees are distinct, these two 
distributions can be combined to produce four different topolo- 
gies, which are defined in Table [T] The probability of any SF 
topology occuring in a randomly constructed network decreases 
rapidly as the number of nodes becomes larger than the average 



connectivity (^Oikonomou and Cluzel 2006 1. 



Since the observation of the Internet's SF topology by ( AI 



bert et al. 1999 1, SF networks have been identified in many 



other fields, ranging from citation networks (^Barabasi and Al- 



bert 1999 1 to biological networks, including metabolic net- 



works Peong et al.|[2000l ) and GRNs fBabu etal .112007 ). Given 
the improbability of large SF topologies occurring randomly, 
this ubiquity has been suggested to be a consequence of com- 
mon generating mechanisms (Barabasi and Albert 1999[ l. One 
possible mechanism is preferential attachment, whereby highly 
connected nodes have a greater probability of receiving a new 
connection — "the rich get richer". This could certainly ex- 
plain the power law distribution of networks such as the Inter- 
net (Barabasi and Albert |1999[ l, but may not be the only re- 
sponsible mechanism ("Keller 2005 1. A modified view of pref- 
erential attachment has been suggested to be responsible for 



the power laws observed in large GRNs ( Barabasi and Oltvai 



pOCMV A fundamental process in biological evolution and GRN 
growth is gene duplication (Carroll 2005 ), where a gene is sim- 
ply copied twice into the daughter system. In the corresponding 
GRN, all connections from the duplicated gene are doubled. As 
a highly connected gene is more likely to be connected to the 
randomly duplicated gene, new genes preferentially attach to 
highly connected genes. 

The process outlined above would give rise to an SF degree 
distribution. However, in GRNs, it is only the outgoing de- 
gree that has been observed to possess a power-law distribu- 
tion (Babu et al. , 2004; jAldana et al.| |2007| ). Either evolution- 
ary rewiring has undone SF incoming degrees (due to the vast 
number of random network configurations relative to scale-free 



ones), or it is not gene duplication that is responsible for the 
generation of power-laws within GRNs. In either case, there 
must be another force either maintaining or generating an out- 
going SF distrbution. 

The effect of the SF topology on Boolean networks has al- 
ready been investigated by other authors. For example, Al- 



dana ( 2003 1 and Aldana and Cluzel ( 2003 1 investigated the ef- 



fect of dynamical perturbations in Boolean networks with an SF 
topology. In random networks, the ordered phase is achievable 
through all nodes having low connectivities, with fine tuning 
of mean degree connectivity. For SF topologies, however, the 
ordered phase can be attained whilst allowing the presence of 
highly connected nodes, or "hubs". In more recent work, Al- 
|dana et aL] ( |2007| ) have also demonstrated increased robustness 
of the attractor landscape in SF topologies under the process of 
gene duplication. In an important paper, [Oikonomou and Cluzel] 
( 2006| l studied large (A^ = 500 node) networks with Boolean 
threshold dynamics (Kurten 1988). They used standard genetic 
algorithms (GAs) to model evolution and found that when the 
selection pressure is for a single node in the network to pro- 
duce a desired oscillatory target, an SFin topology will on av- 
erage reach its target in significantly less generations than an 
ER network will, suggesting that the former is more evolvable 
than the latter. The enhanced evolvability for oscillatory targets 
holds for a number of diff'erent average connectivities. They ar- 
gue that, in contrast to ER networks where the parameters must 
be fine-tuned to achieve a critical phase with its associated en- 
hanced evolvability, the SF network will promote evolution for 
a wide set of different parameter values. 



In this study we consider the work of Oikonomou and Cluzel 



(2006 1 further, with a particular focus on how the enhanced 
evolvability they find depends on the specific target (phenotype) 
and how it relates to both mutational and environmental robust- 
ness. Besides the SFi,, topology used in Oikonomou and Cluzel] 
( 2006| l we also investigate the other topologies in Table [T| with 
a particular focus on the SFout topology that better resembles 
biological GRNs. We note that in our study, as in Oikonomou] 
and Cluzel ( 2006| l, the emphasis is on general results for very 



schematic models of GRNs rather than on networks directly ex- 
tracted from particular experimental systems. The advantage of 
this approach is that general trends are easier to extract, while 
the disadvantage is of course that any direct connections to bi- 
ological systems will need to be worked out further down the 
line. 

We proceed as follows: In Section l2] we describe in some 
detail the methods that we use to simulate the evolutionary dy- 
namics of the networks. In Section[3]we study the evolvability 
of the four different topologies described in Table [T] for single 
oscillatory targets of length L - 10, confirming that all types of 
SF topology are more evolvable than an ER topology. We also 
demonstrate that SF topologies are much more likely to gener- 
ate oscillatory outputs even without any evolutionary pressure. 
Section|4]introduces two new kinds of targets. For multiple os- 
cillatory targets of the same length L, the difference between ER 
and SFout topologies is enhanced for increasing number of in- 
dependent targets. However, these differences are much smaller 
for multiple targets of different length. The difference in evolv- 



ability between the ER and the SFout topologies is also quite 
small when evolving towards a target that fixes 100 of the 500 
output nodes, suggesting that the enhanced evolvability of the 
SFout topology depends in part on the kind of phenotype that 
is being selected for. In Section |5] the mutational and envi- 
ronmental robustness of the ER and SFqu, topologies are stud- 
ied. Even though the SFout topologies are much more evolv- 
able than the ER topologies for oscillatory outputs, they are 
simultaneously more mutationally and environmentally robust. 
Under neutral evolution both the mutational robustness and the 
environmental robustness increase. We discuss our main re- 
sults in Section |6] To show that these results are not limited to 
the model used in Oikonomou and Cluzel| ( ;2006j , we demon- 
strate in Appendix A that another GRN model ( Wagner 1994[ l 
also exhibits a similar interplay between robustness, evolvabil- 
ity and topology. Appendices B - E are shorter and mainly 
focus on technical points left out of the main paper. 

2. Methods 

We model the evolutionary dynamics of a population of 
GRNs in several stages. Firstly, sample networks are created, 
generating a population of genotypes. This is followed by the 
simulation of network dynamics determining the attractor, with 
the phenotype obtained from measurements at a node or set of 
nodes randomly selected at the beginning of the run. This pop- 
ulation of networks then undergoes evolution through recurrent 
mutation of the genotype, followed by selection for the fittest 
individuals based on their phenotype, which is re-measured ev- 
ery generation. 

2.1. Structure 

Networks were created using the "Configuration model", 
a standard method of generating arbitrary directed net- 
works ( Newman et al.| [200T] |Milo et al.| [2504) 1. The num- 
ber of nodes in the networks, A^, was set to 500, whilst the 
average connectivity of network topologies was chosen to be 
1.88. These values provided a direct comparison with the work 
of |Oikonomou and Cluzel|p^006| l, as well as being biologically 
reasonable ( Alon 2006'). 



To allocate node degrees for SF topologies, numbers were 
drawn randomly from a power-law distribution, with the prob- 
ability of degree k given by 



PsFik) 



k-y 



m-y 



1 <k<N 



(1) 



where y parameterises the distribution. 

For the ER topology, in the case where the number of nodes is 
much larger than the average connectivity, the binomial distri- 
bution may be approximated as a Poisson distribution ("Newman' 
|et al. 2001) . In contrast to the SF distribution, where all nodes 
have at least one connection (Psf(O) - 0), the Poisson distri- 
bution permits some nodes to have either an in or out degree of 
0, or even both. Such nodes cannot both affect and be aff'ected 
by the network and thus, their presence introduces a reduction 
in the effective size of the networks. To combat this effect, we 



imposed a condition on the Poisson distribution, whereby de- 
grees of are disallowed and the distribution is re-normalised 
accordingly. We therefore restrict ourselves to considering net- 
works consisting of a single, large connected component. This 
adapted Poisson distribution for a degree k is given by 



PER(k) = 



K'e 



k„-K 



(1 -e-'^)k\ 



I <k<N 



(2) 



where K parameterises the distribution. 

The average connectivity of each distribution can be deter- 
mined simply by (k) - Yj^^i kP(k). The corresponding values 
of 7 and K for an average connectivity of 1.88 are y - 2.5 and 
K ^ 1.431. 

2.2. Dynamics 

A Boolean threshold model was used to model the dynamics. 
At the beginning of a dynamical run, each node o", is initialized 
by a random process to to be either in an "on" state (cr, = 1) or 
an "off" state (cr, = 0). These are the initial conditions (ICs) of 
the network. In an evolutionary run, we either impose the same 
set of initial node states at each generation (constant ICs) or 
choose initial node states randomly at each generation (random 
ICs). We investigate both situations in this work (in [Oikonomou 



and Cluzel (20061 only random ICs were used). From an ini- 



tial configuration, the network dynamics are iterated over a set 
of discrete time steps, with the state of each node updated syn- 
chronously between time steps (see Fig. [T^). 

The rules for updating the state of a particular node depend 
on the state of the incoming nodes, and the strength of these 
signals. An NxN matrix of weights, w, defines the interactions 
between all nodes. Connections between nodes are assigned 
weights in matrix w randomly from the set of real numbers on 
the interval (-1, 1). 

The state of node / at the following time step is denoted cr,(f-i- 
1) and is determined by a sum over all incoming nodes 






(3) 



combined with the following threshold rule: 

cr,(f+l) = 



1 if5,>0 

o-,(0 if Si ^0 
ifSi<0 



(4) 



To physically interpret the matrix of weights, consider a con- 
nection from node ; to node j. If w,y > 0, node j promotes node 
; whilst for w,-; < 0, node j inhibits node /. If w/j - there is no 
connection from J to ;. 

To prevent frozen dynamics, there is a modified rule for 



nodes with only a single incoming connection ( Oikonomou and 



Cluzel 2006 1. If node j is incoming to node i the the dynamics 
depend purely on the sign of the weight with 



\o-j(t) ifw^-,>0 
l-icr/f) lfw;,<0. 



(5) 
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Figure 1 : A summary of the dynamical and evolutionary processes acting on a network. In (a) an example of the network dynamics is demonstrated. The node 
coloured blue is selected as the output node. When one time step is performed, the state of the nodes will have been updated synchronously by the rules associated 
with Eqns. [3]and|4] Further times steps are then performed to produce an output bit string as shown above. This is then compared to the target and its fitness 
calculated with Eqn. |6] In this particular case maximal fitness is achieved. In (b) the two types of mutation are demonstrated. The red arrows represent a swapping 
mutation where the grey edges have been switched but preserve each node's degree distribution. The green edge, with a new weight of 0.79, represents a weight 
mutation from its previous value of 0.23. When a mutation occurs during an evolutionary run, one mutation type is chosen, with each type being equally probable. 



where -> is the NOT operator. The behaviour of nodes con- 
nected by this modified rule is simple. Node states are either 
copied or inverted for positive and negative weights respec- 
tively. 

The critical parameter values, jc and Kc, separating the or- 
dered and chaotic phases for each topology, have been previ- 
ously calculated by Oikonomou and Cluzel| (|2006 ) based upon 



2.3. Evolution 



To model evolution of the population towards a desired out- 
put (phenotype) we used a simple genetic algorithm. First a 
population of Npop - 50 independent networks is generated. 
Then, for each generation, the following cutoff selection regime 
is applied. Each network first produces three daughter net- 



the ann ealed approximation introduce d by Derrida a ndPomeauj ^Qj-ks, enlarging the population from Npop = 50 to A^^^^, = 200. 



( 1986 1. We repeat the calculation (see Appendix B i finding the 
critical value for the SF topology to be jc - 2.42 for networks 
of size A^ - 500. We also calculate the critical connectivity 
for the adapted Poisson distribution, which we determine to be 
Kf- - 3.538. The primary parameter values used in this study 
of 7 = 2.5 and A" = 1.431, both produce average connectivities 
below that needed for criticality, so that the networks are within 
the ordered regime. 



After a run typically lasting 350 time steps, the steady state 
dynamics of the output node(s) in each network is determined 
and this sequence defines the individual's phenotype. An ex- 
ample of a single dynamical time step is performed on model 
networks in Fig. [Tk. 



Each node in these daughter networks has a mutation proba- 
bility fi = 0.02 of undergoing a mutation, a value previously 



used by Oikonomou and Cluzel (2006 1. There are two types 
of mutation: either a random change in the weight of any of a 
node's connections or an incoming connection is swapped with 
another incoming connection in the network, provided the new 
connections do not already exist (Fig. [T]i. Both mutation types 
preserve the node degree distribution throughout an evolution- 
ary run. 

In single target experiments, the output of a node, randomly 
selected at the start of the evolutionary run, is used to calculate 
the network's fitness. The fitness is calculated from the min- 
imum Hamming distance between the node's periodic output 
and the target output over all cyclic permutations. We define 



this mathematically as 



F = max 

je|0,...,<5) 



1 -''-KK ' 



'-'o'^tgt "J^ 



(6) 



where Lo is the output period, L,g, is the target period, cro(i) is 
the ;th bit in the output sequence and (Ttgtii) is the /th bit in the 
target sequence. The sequences are, for comparison, extended 
by repeating them periodically to construct sequences of equal 
length LoLtgt- The sum is then a Hamming distance of these two 
strings. Taking the maximum over an offset j between and 
S - max{Lo,L,s,,) performs the Hamming distance calculation 
over all cyclic permutations of the strings. 

The Hamming distance measures the fraction of the sequence 
that matches the target, providing a linear measure for the fit- 
ness. Once all networks in the enlarged population have been 
assigned a fitness, the networks are ordered by fitness and the 
top 50 kept to form the population for the next generation. 

3. Evolution towards single L=10 targets 

3.1. All types ofSF topology are more evolvable 

Expanding on the work of |Oikonomou and Cluzel|f2006| l, we 
tested evolutionary performance or evolvability of each of the 
four different topologies (ER, SFout, SFi„ and SF). Evolvability 
is measured through adaptation speed towards the target output, 
i.e. the number of generations before every member of a pop- 
ulation has reached maximum fitness. Evolvability has been 
frequently used in different contexts, sometimes creating con- 



fusion ( Poole et al. 2003 Wagner 2005) . In this case, however, 
the adaptation speed would seem the most natural measure of 
evolvability. 

A randomly chosen node from each network was evolved 
towards a randomly chosen L = 10 target for a maximum of 
10000 generations. This was repeated 50 times for each topol- 
ogy, with both constant and random ICs. When evolution is per- 
formed with random ICs, robustness to this stochasticity must 
evolve for an individual to maintain itself in the population. For 
constant ICs there is no such constraint. 

We observed that occasionally an evolutionary run continues 
for a very large number of generations before converging to a 
solution, or sometimes not converging at all within the limits of 
our calculations. The mean of the adaptation time distribution 
can be dominated by these rare events. For that reason, we use 
the median adaptation time, arguing that this provides a more 
accurate measure of the typical adaptation time than the mean 
does. We define t to be the number of generations for all mem- 
bers of a population to reach maximum fitness in a single evo- 
lutionary run, whilst f is the median over a batch of runs. The 
values of f for the different topologies and ICs, are presented in 
Table 121 

Out of the three SF topologies, the SFboth evolves more 
rapidly than the SFin topology does, which in turn is fol- 
lowed by SFout and ER. This result supports the findings 
of|Oikonomou and Cluzel| ([2006[l, where the SFi,, topology is 



Table 2: The median number of generations to reach maximal fitness, f , for runs 
evolving towards a single L = 10 target. Results are shown for both random 
and constant ICs. All SF topologies outperform the ER topology. The SFjn 
topology is most significantly aft'ected by the dift'erence between constant and 
random ICs, whilst the SFout is least affected. 
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SFout 


SFu, 


SF 
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3 200 
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tested against the ER topology and shown to evolve much more 
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Figure 2: The probability of different periods being the largest in randomly 
generated networks. SFboth and SFj^ perform similarly, whilst SFout and ER are 
similar too. 



rapidly. The order of the topologies, with respect to adapta- 
tion speed, is unchanged between constant and random ICs, 
although the decrease in speed from constant to random ICs 
varies amongst the topologies, with SFin and SFboth topologies 
showing greater relative changes than ER and SFout do. This 
latter property gives an indication of the natural robustness to 
ICs of each topology, a property we explore in later sections. 



3.2. The SF topology exhibits more oscillations and greater 
synchronicity than the ER topology 

To investigate whether the SF topology more naturally gen- 
erates oscillations than the ER topology does, we examined the 
dynamics of networks without evolution towards a target. For 
each topology 10000 networks were generated and two mea- 
surements made. For the first, the period of every node in a net- 
work was found and the largest period recorded. A histogram 
of these periods is presented in Fig. |2] For the second measure- 
ment, the period of a randomly chosen individual nodes was 
recorded. A histogram of periods from these 10000 networks 
is presented in Fig. [3] 

In Fig. [2] we observe a clear divide between the SFout topol- 
ogy and the other two SF topologies. The SFout topology pro- 
duces a probability distribution similar to that of the ER topol- 
ogy. The SFboth and SFju topologies have similar distributions 
to each other They are significantly more likely to have maxi- 
mum periods at larger values of L. 
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Figure 3 : The probability of different periods being measured at a random node, 
in randomly generated networks. Note the logscale. SFboth generates the most 
oscillations, followed by SFin, SFout and ER. 



Fig. l3]examines the probability to find oscillations of period 
L at a randomly chosen node. This result differs to the previous 
one as SFboth now shows more oscillations than SFin, whilst 
SFout is clearly superior to ER at generating oscillations. Again 
the SFjn and SFboth topologies are considerably more likely than 
the SFout topology to have an oscillating node. In both Fig. l2] 
andlSleven oscillations are more easy to generate than odd ones. 

We next measured the distribution of oscillations in the net- 
work during evolutionary runs towards L - 10 targets. The 
period of each node was measured in the first generation where 
all individuals in the population have a fitness of unity (defined 
to be the first maximally fit generation), giving an indication of 
the extent to which the target node's period is spread throughout 
the network. This measures the synchronicity within evolved 
networks. 

The mean fraction of nodes oscillating with each period are 
presented in Fig. HI for each of the four network topologies, 
averaged over 50 independent evolutionary runs. This figure 
demonstrates that a greater fraction of ER networks are frozen 
{L - 1). All the SF topologies have an increased number of os- 
cillatory nodes and particularly L = 10 nodes. Periods that are 
factors of 10 (2 and 5) also appear to be much more prevalent 
than other periods. 

In general we observe that the scale-free topologies are more 
likely to exhibit oscillations at a randomly chosen node, have 
larger maximum periods and show more synchronicity during 
an evolutionary run. In general the differences with the ER net- 
works are largest for the SFi,, and SFboth topologies. In that con- 
text it is important to recall that there is a modified update rule 
for single incoming nodes. Through this rule, any oscillatory 
node connected to a node with a single incoming connection 
will always propagate the oscillation to that node. Given that 
Psf(_1) - 0.745 compared with ^£^(1) - 0.450, guaranteed 
propagators are much more likely in SFin networks. The adap- 
tive advantage generated by the SFout topology must have dif- 
ferent origins, however, as the corresponding P{\) is the same 
as for ER networks, and hence the modified rule affects both 
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Figure 4: The probability of a randomly chosen node possessing period L in 
the first generation where maximum fitness for an L = 10 target is attained. 
The probability is calculated by taking the mean fraction of nodes with each 
period, taken over 50 independent evolutionary runs. The "oth" column refers 
to other periods, i.e. all periods not labelled. This figure measures the degree 
of synchronicity in the networks. Note the difference in scale on the y-axes of 
the right and left graphs. 



cases similarly. 

In the rest of this paper we focus mainly on comparisons 
of the ER with the SFout distribution because the latter most 
closely resembles the topology of GRNs found in nature (Babul 
|etaLl|2004l[Aldana et al.||2007] l. 



4. Evolution towards multiple oscillatory and frozen tar- 
gets 

In nature GRNs may employ phenotypes where multiple 
nodes require oscillating periods. This would be the case, for 
example, if a set of genes played a part in an oscillatory pro- 
cess within an organism. The segmentation clock (Pourquie 



2003) and circadian rhythms ( |Hirata et al.||2002} rely on multi- 
ple genes working together in this way. 

To study this scenario, evolution was performed towards a 
target phenotype defined by requiring the output of several ran- 
domly chosen nodes' to have specific periods. We define the 
target period set, £., to comprise the period of each target node. 
For example, for A^tgts = 3 with two L - 5 targets and one 
L= 10 target, X = {5,5, 10). 

To determine the fitness of a network, each target is assumed 
to make an equal contribution. A random set of nodes are as- 
sessed against each of the randomly chosen target outputs of 
length specified in £. and each node's fractional fitness is added. 
As such, we define fitness of the networks with multiple targets 
as 



A'lsi 



2^'/^'g'^ 



(7) 



/=i 



where Fi is the fitness of the /* node in the phenotype and A^tgts 
is the number of target nodes defining the phenotype. Each F, 
is given by Eq. (|6]l. 
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Figure 5: f v.s. the number of L = 7 targets for the ER (red, X markers) and 
SFout (blue, + markers) topologies with corresponding fits (the lines are guides 
to the eye). Both topologies experience an increase in adaptation time with 
an increase in the number of targets. This increase is more rapid for the ER 
topology, typically requiring more than 20 000 generations for Ni^ts > 6. The 
slower rate of increase of f seen in the SFom topology is likely to be due to the 
greater synchronicity of this topology. 



4.1. The adaptation time required to achieve multiple targets 
of the same length increases rapidly with the number of 
targets 

As the ER topology evolves rather slowly towards L = 10 tar- 
gets, a smaller target length was selected. This allows a greater 
range of target numbers to be tested within a computationally 
reasonable number of generations. If the target length is too 
small, however, only a few different targets will be possible. 
For example, there are only 6 possible L - 5 targets (see Ap 
Ipendix C"| . Given these considerations, multiple L 



1 targets 

were chosen as this length is easily attainable by the ER topol- 
ogy whilst there are 18 possible L-1 targets. 

To measure speed of adaptation, we used the median number 
of generations until a population is maximally fit, f, as in Sec- 
tion 3.1 In Fig. [5] T is plotted against A^tgts- As expected, for 
both ER and SFout topologies there is an increase in the adap- 
tation time as A^tgts increases. The change in adaptation time is 
non-linear, increasing rapidly with the number of targets. The 
rate of increase is much more rapid in ER networks than SFout 
ones. For example, SFout networks are still capable of adapting 



to Mgts 
A'tgts = ^ 



20 in just over half the time ER takes to adapt for 



As observed earlier in Section 3.2 the SFout topology shows 



greater synchronicity for oscillatory signals than the ER topol- 
ogy does. This difference may explain why the SFout topology 
performs so much better for multiple targets of the same length. 

4.2. The adaptation time required to acheive multiple targets 
of different lengths depends on the number of targets and 
their relative periods 

To test whether the enhanced synchronicity of the SFout 
topology is a key reason it performs so much better for multiple 
targets of the same length, we also examined different sets of 
multiple targets with different lengths. All the evolutionary runs 



Table 3: Values of f for evolutionary nins evolving under random ICs to differ- 
ent sets of periods, f is smaller for £. = {5, 101 compared with £. = {5, 8} for 
both topologies. For the X = {5, 10, 20) no convergence was achieved. 



Period set, £. 


ER 


SFout 


{5,10) 

{5,8) 

{5,10,20) 


7 000 
9 900 
> 40 000 


2400 
4700 
> 40 000 



were performed with random ICs, and the results are shown in 
Table l3] Clearly having targets with multiple lengths greatly 
slows down the median adaptation time compared to having 
multiple targets of the same length for both topologies. It is 
also interesting to compare the performance of the networks for 
the £. - {5, 10) targets to the performance for the £. - {5, 8) 
targets. Although the former targets are slightly longer than the 
latter, they are easier to reproduce. Presumably this is because 
L - 10 is a multiple of L = 5 whereas L = 8 and L - 5 are co- 
prime. We also note that the difference between the SFout and 
ER topologies is smaller for the coprime targets, suggesting that 
the enhanced synchronicity of the SFout topology, observed in 
Fig. HI plays an important role in its enhanced ability to adapt to 
certain types of targets (phenotypes). Finally, we also tested the 
target X = {5, 10, 20) but were not able to find solutions within 
the 40 000 generation cutoff we used in our simulations, show- 
ing that multiple longer targets are significantly more difficult 
to achieve than a larger set of shorter targets (as in Fig.lSll. 



4.3. The evolv abilities of SFom and ER networks for multiple 
L — I targets are much more similar than for oscillatory 
targets 

Thus far only targets of an oscillatory nature have been con- 
sidered. This is because single L-1 targets would be trivial 
to evolve; nodes of L = 1 are by far the most common in all 
attractors. With the extension to multiple target nodes, how- 
ever, evolving to a large number of specific L - I targets is 
a non-trivial task. Biologically, a set of "frozen" targets could 
correspond to a particular constant gene expression pattern over 
some length of time, for example in an organisms' develop- 
ment. 

As was done earlier, networks of each topology were evolved 
over 10000 generations in 50 independent runs. Out of the 
A^ = 500 nodes, 100 were randomly chosen and required to 
possess a specific L = 1 target. We will refer to this target set as 
X. - {100 X 1). The values off with random and constant ICs 
are displayed in Table HI The SFout topology maintains its ad- 
vantage over the ER topology but the difference in evolutionary 
adaptation times are vastly reduced. The effect of random ICs 
does not slow evolution down that much, and there is a similar 
relative decrease for both topologies. 

5. Robustness in single & multiple target experiments 

In the previous section we have investigated the evolvability 
of ER and various SF topologies towards oscillatory and frozen 



Table 4: Median adaptation time f for runs evolving towards X = {100 X 1| 
targets. The performance of the topologies is much more similar than for evo- 
lution towards oscillatory targets, although the SFout topology retains a slight 
advanatage. 



ness, R„,. By averaging over the population, the statistical error 



is reduced. We show in Appendix E that the diversity of geno- 
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Figure 6: Mean mutational robustness R„,, against generation T (where T = 1 
is the first maximally fit generation) for X = {10|. The errorbars represent the 
standard error on the mean. Robustness is still increasing after 10000 gener- 
ations for both topologies. The SFom topology is more robust from the outset 
and remains so throughout. Both topologies are able to evolve greater robust- 
ness. Noting the v-axis range, the change in robustness is small but statistically 
significant for both topologies. 



targets. In this section we study the mutational and environ- 
mental robustness of these networks under similar evolutionary 
situations. We particularly focus on the change in robustness 
under neutral evolution, after the population has achieved max- 
imal fitness. 

5.7. Mutational robustness is higher in SF topologies and in- 
creases under neutral mutation 

To examine the mutational robustness of large networks 
with a large range of weights, a sampling approach must be 
taken as the number of possible network genotypes in a 1- 
mutation neighbourhood is astronomical. To sample these 
genotypes, each node in a network was tested with 6 sam- 
ple mutations of each mutation type in Fig. [T] We em- 
ploy constant ICs throughout, ensuring that IC robustness ef- 
fects do not play a role. An individual's mutational robust- 
ness, R'^^, is defined as the proportion of sample mutations 
that do not result in a network's phenotype changing. The 
robustness testing was performed on every member of a pop- 
ulation during an evolutionary run at the generation where 
maximum fitness is attained (T - 1), and then at intervals 
(T = 2,5,10,20,50,100,200,500,1000,2 000,5 000,10000) 
afterwards. At each interval, R"f is averaged over the whole 
population producing an average population mutational robust- 



types is maintained during neutral evolution. 

Mutational robustness was measured for X = {10) and X. = 
{100 X 1) targets, although we only show results for the former 
target in Fig. [6] The main and perhaps unexpected result for 
both types of target is that throughout these evolutionary peri- 
ods the SFout topology maintains a greater robustness to ran- 
dom mutation than the ER topology, even though it is simul- 
taneously more evolvable. We also note that R,,, is fairly high 
for both topologies, suggesting that the vast majority of mu- 
tations don't affect the phenotype. The most likely reason for 
this is that a large number of nodes simply do not participate in 
the core network that drives the output of the target node. Al- 
though it is not shown, the value of/?,,, is also fairly high for the 
X. = {100 X 1) targets {R,„ > 0.8) and moreover the SFout topol- 
ogy outperforms the ER topology here as well. Again, since 
20% of the nodes are fixed by the target output, this suggests 
that a large number of the nodes do not participate in the output 
that determines the phenotype of these systems. 

The second effect this figure clearly shows is that the mu- 
tational robustness increases with evolutionary time. It was 
shown in van Nimwegen et al. ( 1999| l that under neutral mu- 
tation, the mutational robustness of a population can increase. 
One mechanism for this phenomenon is that those genomes that 
are more mutationally robust are more likely to have offspring 
that share the same phenotype. Under continued selection then, 
such genomes are more likely to survive than genomes with low 
mutational robustness. These authors also showed that for se- 
lection to optimize mutational robustness in neutrally evolving 
systems, the condition fiNpop s> 1 must be satisfied. Here the 
mutation rate per node is set to fi„ode = 0.02. For an A^ = 500 
node network, this gives an average of fj. - ^i„odeN - 10 
mutations per mutated offspring. With a fraction m - 0.75 
mutated offspring in the extended population before the cut 
off is imposed, in each generation there are fi x m x Np„p = 



pop 



10 X 0.75 X 50 = 375 mutations in a future gen- 
eration's genepool. Given that the condition jiNpop » 1 is 
satisfied, along with the genetic diversity observed in the pop- 



ulations (see Appendix E i, further optimisation of mutational 
robustness is certainly expected. However, selection may not 
be the only cause of the observed effects. It may be possible for 
the increase to be obtained purely through entropic effects - for 
example, if many more genotypes are available at higher muta- 
tional robustness, the population will discover such possibilities 
simply through diffusing over the neutral space. Working out 
the exact reasons of this optimisation is thus a subtle question, 
and we don't think that our results unambiguously fix the cause 
of the increase in /?,„. 

5.2. Environmental robustness increases markedly with time 
after maximal fitness is reached for random initial con- 
ditions, and is higher in SF topologies 
GRN oscillators in vivo are often required to function un- 
der a variety of environmental conditions. By possessing en- 
vironmental robustness, function can be maintained under such 
stochasticity and is thus a desirable trait. Here we examine the 




10,000 




10,000 



Figure 7: Mean environmental robustness Rj, against T witli random ICs for 
£. = {10|. SFout is initially much more robust, although the gap is reduced 
over the course of neutral evolution. The increase in Rj is dramatic for both 
topologies, indicating a strong optimisation. 



change in environmental robustness over time during neutral 
evolution and consider runs with both random and constant ICs. 

Given the size of the networks (there are 2^"" ^ 3 x 10'^'^ 
different initial configurations) a sampling approach was nec- 
essary. To estimate the environmental robustness, a random 
set of initial node states was taken, and the output phenotype 
was measured. This process was repeated for 1 000 different 
initial conditions (chosen randomly for each dynamics simu- 
lation) and the average fraction of initial conditions where the 
target phenotype was produced for an individual is called R"''. 
The average of /?™'' was taken over all the networks in the pop- 
ulation to give Rd, a measure of the environmental robustness 
(the subscript d refers to the dynamic nature of this measure of 
environmental robustness). This quantity was then further av- 
eraged over 50 evolutionary runs to reduce statistical error, and 
measured at the same intervals as in Section lSTI after maximum 
fitness is attained. 

Results for the mean environmental robustness R^ for X. = 
{10) targets are shown for random ICs in Fig. IT] and for runs 
with constant ICs in Fig. IS] Firstly we observe that the SFout 
topologies are more environmentally robust than the ER topolo- 
gies are, an effect that mirrors the mutational robustness. Sec- 
ondly, the runs with random ICs show a marked increase of Rj 
during neutral evolution. The runs with fixed ICs show a mod- 
est increase of Rj during neutral evolution, especially for the 
ER networks. 

When random ICs are used, it may be possible for selec- 
tion to increase the environmental robustness of the population, 
for similar reasons to the arguments for the increase in muta- 
tional robustness. During the evolutionary runs those genotypes 
that are more environmentally robust are more likely to have 
progeny that generate the required phenotype. Thus they are 
more hkely to survive than less envrionmentally robust geno- 
types. Furthermore, we can make an estimate of the minimal 
environmental robustness (/?'"'") for individuals to simply be 
sustained in the population. Let us assume individuals possess 
maximal mutational robustness and mutations affect their en- 



Figure 8: Mean environmental robustness R^i versus T, with constant ICs for 
X = 1 10). The SFout topology naturally possesses greater environmental ro- 
bustness. Over time the robustness of SFom changes insignificantly. The ER 
topology, on the other hand, almost doubles its robustness over the period mea- 
sured. However, the ER topology still remains much less robust than the SFom 
topology. 



vironmental robustness negligibly. Under these conditions, at 
least one of the four genotypes must produce the target pheno- 
type again. This gives a value of R""" - 0.25. As the solution 
must propagate through the population, and is not fully muta- 
tionally robust, this is the lower bound and the actual value of 
Rd required may be significantly larger. 

For constant ICs there is no minimum environmental robust- 
ness which explains in part why the values of Rd are lower than 
for random ICs. Moreover, there is no direct selection pres- 
sure for environmental robustness that acts on runs with con- 
stant ICs. We observe that Rd does increase for ER networks. 
The reasons for this are not clear. A possible cause could be 
a correlation between mutational robustness, which is affected 
by selection, and environmental robustness. Another possibility 
would be an entropic effect where diffusion on a neutral space 
is more likely to find phenotypes with higher robustness. 

Although they are not shown, the X = { 100 x 1 ) targets show 
similar trends to the X = 10 targets: the SFout topology is more 
environmentally robust than the ER topology and furthermore 
both topologies show a modest increase in Rd under neutral mu- 
tation for both random and fixed ICs. 



6. Discussion 

In this paper we have studied the interplay between evolv- 
ability and robustness of Boolean threshold dynamics models 
with ER and SF topologies. Our main conclusions are: 

1) The enhanced evolvability of SF networks over ER de- 
pends on the phenotype under selection. 
We confirm the results of [Oikonomou and Cluzel| ( |2006[ ) that 
networks with SFin topology are more evolvable towards sin- 
gle oscillatory targets than ER networks are, and find the same 
enhanced evolvability for the SFboth topology and also for the 
SFout topology that most closely resembles the global topology 
of biological GRNs. The evolvability advantage of the SFout 



10 



topology over the ER topology increases significantly for mul- 
tiple targets with the same length. By contrast, for multiple tar- 
gets of a different length, the difference in evolvability between 
SFout and ER topologies is reduced, and if the targets are co- 
prime, the differences are even smaller. We also studied targets 
where 100 of the 500 nodes were given a fixed output. Here 
again the differences in evolvabilities between the topologies 
was greatly reduced compared to the L = 10 oscillatory outputs 
studied in |Oikonomou and Cluzel| ( |2006| l. These results show 
that the evolvability advantage of a particular network topology 
depends strongly on the phenotype under selection. 

2) SF networks show more oscillatory phenotypes and 
greater synchronicity than ER networks do. 
We compared the likelihood of finding an oscillatory output of 
period L i^ 1 at a randomly chosen output node in networks 
without any evolutionary selection. The SF networks show a 
greater probabiUty of finding oscillatory outputs than the ER 
networks do. Similarly, once a population has reached max- 
imum fitness under evolution towards an L = 10 target for a 
single output node, a different node in the SF networks is much 
more likely than an ER network to show L = 10 oscillatory out- 
put: SF networks show much greater synchronicity than ER net- 
works do. We argue that the increased evolvability towards os- 
cillatory targets observed for SF networks is due to the fact that 
such phenotypes are more naturally generated by these topolo- 
gies. 

3} SFout networks are more mutationally robust than ER net- 
works are. For both network topologies, the mutational ro- 
bustness increases under neutral mutations. One might naively 
expect that the enhanced evolvability would lead to a reduced 
mutational robustness. However, as has been pointed out by 
numerous authors ( |Bloom et al.| |2007a|b[ [Aldana et al.| |2007[ 



We also find that the mutational robustness increases un- 
der neutral mutations, an effect we attribute to the fact that 
more mutationally robust genotypes are more likely to have 



Wagnerl |2008b"a; 'Daniels eFaLl [2008} [Draghi et al.|[20T0^, this 



simple expectation does not necessarily hold for populations. 
Even without selection, SF topologies show a greater propen- 
sity than ER networks to exhibit oscillatory phenotypes. This 
difference suggests that the number of different genotypes that 
generate the same oscillatory phenotype is larger in SF net- 
works, i.e. that the neutral spaces are larger Larger neutral 
spaces may explain how greater mutational robustness can cor- 



relate with greater evolvability (Wagner 2008b I. However, to 
firmly establish that the evolvability advantage of the SF net- 
works can be explained by properties of the neutral space re- 
quires further investigation. For example, one must show that 
larger neutral spaces have evolutionary access to a greater phe- 
notypic diversity. 

Another possible factor is the heterogeneity introduced by 
the SFout topology. The measure of mutational robustness treats 
each node in the network equally (a biologically reasonable as- 
sumption), whereas, in an SFout topology, nodes are not equal: 
most nodes trivially have a single outgoing connection, whilst 
a few have a large number, the so-called 'hubs' of the network. 
The hubs in the SFout topology may play an important role, and 
mutations to these may be most effective, but rare as hubs only 
comprise a small part of the network. However, when muta- 
tions do occur to hubs (or to nodes closely connected to a hub), 
there is the potential for large scale influence of the mutation. 



progeny that survive under future mutations ( .van Nimwegen 
|etaLi[T999l l. 

4} SFout networks are more environmentally robust to 
changes in initial conditions than ER networks are, and both 
networks types exhibit an increase in environmental robustness 
under neutral mutations. The increased environmental robust- 
ness observed in SF topologies may be due to the comparatively 
large number of nodes acting in cohort to produce the required 
oscillatory signal - a feature shown by the observed high syn- 
chronicity in SF topologies. This large group of synchronised 
nodes may act to increase the number of initial states that con- 
verge to the required attractor, as perturbations to node states 
may be "damped out" by the overall tendency of the group to 
produce a particular oscillation. By contrast, the smaller groups 
of synchronised nodes in ER networks may be more sensitive 
to the initial states of nodes. 

For random ICs, the environmental robustness increases un- 
der neutral mutations for reasons similar to those that explain 
how mutational robustness increases: more environmentally ro- 
bust genotypes are more likely to survive in future generations 
when the environment changes. We find a similar, but more 
modest increase in environmental robustness in ER networks 
for fixed ICs, where there should be no selection for more ro- 
bust phenotypes. At present it is not clear what causes this in- 
crease, but it may be correlated to the increase in mutational 
robustness. 

An important question is whether the conclusions listed 
above hold only for the particular models we studied, or 



whether they are valid more generally. We show in Appendix 



lAI that similar conclusions also hold for a different threshold 
model ( Wagner||199 4) that is frequently applied to study GRNs. 
Nevertheless, the evolutionary dynamics of the threshold mod- 
els we study here is considerably simpler than that which oc- 
curs in nature and so further work is needed to probe to what 
extent these conclusions hold for biological GRNs. At present 
such studies are very challenging because they require knowl- 
edge of the evolutionary history of GRNs, much of which is 
not easily accessible. The conclusions from this work that will 
most likely carry over are: (a) that the relationship between 
robustness and evolvability can depend on the class of pheno- 
types under selection; and (b) an organism's GRN can become 
both more evolvable and more robust simply though changing 
its overall topology. 
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Appendix A. Comparison with the Wagner model 



To examine the generality of the dynamical properties intro- 
duced by the SFout topology, we made a comparison with an- 
other threshold model that was proposed by [Wagner ( 1994[ l and 
that has been investigated extensively with respect to GRNs. 
The essential differences between the Wagner model and the 
one of [Oikonomou and Cluzel] ( |2006| l, is the use of different 
node states and update rules. In the original Wagner model, the 
gene states for node / are either ctj - 1, 0, -1. The update rules 
differ by using purely threshold rules, with no modified rule for 
single incoming degrees. The node state ctj - 0, occurs when 
the input sum S t - 0. Node states of 1 and - 1 correspond to our 
"on" and "off" states considered in the lOikonomou and Cluzell 
(|2006 ) model. 



We investigated dynamics with the essence of the Wagner 
model. We note that nodes with o-j - make no contribution to 
the dynamics in the next time step and given a continuous range 
of weights, such nodes will only occur when lacking any in- 
coming connections. As a result, nodes with ctj - will always 
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remain in this null state throughout a dynamical run, render- 
ing them ineffective within a network. Effectively, the update 



rules only differ from the Oikonomou and Cluzel (2006) model 



because we can abandon the modified rule for single incoming 
connections. The effect of this rule is now incorporated into the 
behaviour of all nodes, as when a node is in its "off" state, a 
signal is still affective on the nodes it connects to. This is due 
to the "off" state now corresponding to cr,- = -1, resulting in 
a non-zero amount being added in the threshold sum S ,, even 
when the node is "off". 

To examine the networks using the Wagner model, we firstly 
consider for what parameters the system enters the chaotic 



phase. In the Oikonomou and Cluzel ( 2006 j ) model, networks 
were within the ordered phase and need the same to be true here 
for comparison and biological realism. It was found that for the 
adapted ER topology (or SFin topology), a phase transition oc- 
curs for an average connectivity of (k) = 1 . This gives a critical 
parameter of Kc — > oo for the adapted Poisson distribution. As 
a consequence, networks with this topology (or SFin) cannot be 
produced in the ordered phase. This is an effect due to local 
noise propagation in the Wagner model being more likely than 



in the [Oikonomou and Cluzel, (2006^ model (see Appendix B i. 
As a result we considered ER graphs with a standard Poisson 
distribution 

PERjk)^^^ (A.l) 

rather than an adapted one. For the standard Poisson distribu- 
tion, the phase transition occurs at Kc - 2.075. This allowed 
us to compare non-adapted ER networks and SFout (SFout also 
possesses a non-adapted ER incoming degree distribution here) 
networks at the previous average degree connectivity, whilst re- 
maining in the ordered phase. 

Appendix A.l. The evolvability of SF networks using the Wag- 
ner model is again greater than that of ER net- 
works 



As in Section |3.1| we evolved networks towards random 
L = 10 targets over 50 evolutionary runs. The median gen- 
eration numbers for a population to reach maximum fitness, f. 



are shown in Table A. 5 for each topology under random and 



constant ICs. Evolution with constant ICs is a much simpler 
task for both, with t reduced by more than an order of mag- 
nitude for both topologies. Interestingly, the SFout topological 
advantage has vanished in this case. For random ICs, however, 
the results are dramatically different. The SFout topology now 
outperforms the ER topology by roughly the same factor as 



in the Oikonomou and Cluzel (2006) model. Both topologies 
struggle to a much greater extent to adapt under random ICs. 
This is likely to be a consequence of local noise propagation 
being greatly increased in this model. 

Appendix A.2. Networks with the Wagner model show more os- 
cillatory behaviour 



As in Section 3.2 the synchronicity in networks under evo- 
lution towards a single L = 10 target was examined but with the 
Wagner model. The period of each node in the first maximally 
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Figure A. 9: The probability a random node has a given period, L, in the first 
generation where maximum fitness is attained, within the Wagner model. The 
probability is derived from the mean fraction of measured nodes taken over 50 
independent evolutionary runs. The "oth" column refers to other periods, i.e. 
all periods not labelled. As in the model of |Oikonomou and CluzelH2006) the 
SFout topology exhibits greater synchronisity. 



fit population was measured, with an average taken over 50 evo- 
lutionary runs. A comparison of the probability of measuring 



each period is presented in Fig. A. 9 As for the Oikonomou and 



Cluzel ( 2006 1 model, the probability of a node having the same 
period as the evolved output (L = 10) is significantly higher 
in the SFout topology than in the ER topology. In the Wagner 
model, oscillations spread much more widely with both topolo- 
gies having a much smaller fraction of L = 1 nodes. This is 
due to an increase in local noise propagation under the Wagner 
update rules. 

Appendix A.3. Mutational and environmental robustness of 
networks with the Wagner model also increase 
with time after maximal fitness is reached 

In Section |5] we considered the mutational and environmen- 
tal robustness of networks evolving towards £. - {10) tar- 
gets. These experiments made use of the update rules used 



by |Oikonomou and Cluzel ( 2006) . In this section, we confirm 
the overall results presented previously are not a consequence 
of those update rules, and are a general feature of threshold 
networks. The same method of measurements were used as in 
Section|5] 

Fig. |A.10| shows the average change in mutational robustness 
after the population reaches maximal fitness. We find, once 
again, that SFout topologies exhibit greater robustness through- 
out and that both topologies increase their mutational robust- 
ness over time. For the Wagner model the mutational robustness 
changes to a greater extent for both topologies than is observed 
for the model of |Oikonomou and Cluzel |f2006[ ). This difference 
is likely due to the greater extent of local noise propagation in 
the Wagner model. 



Fig. A. 1 1 and Fig. A. 12 show the change in average environ- 
mental robustness over the maximally fit evolutionary period, 
for both random and constant ICs. As before, in the case of ran- 
dom ICs, environmental robustness can be optimised directly 
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Figure A. 1 0: Mean mutational robustness R^ , against T for evolution of L = 10 
targets in the Wagner threshold model. As found for the |Oikonomou and Cluzel| 
^006^ model, neutral evolution can optimise muational robustness over time. 
The relative change in mutational robustness attained in this case is larger than 
that in the |Oikonomou and Cluzel| j20Q6i) model, whilst the typical maximal 
mutational robustness values attained after 10000 generations of neutral evolu- 
tion are similar. 
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Figure A. 11 : Mean environmental robustness Rj, against T for evolution of L = 
10 targets in the Wagner threshold model for random ICs. Selection increases 
environmental robustness in both topologies. The limit of this increase appears 
to be reached towards the end of the neutral evolutionary period, levelling off 
at around Rj » 0.6. The overall lower level of environmental robustness is a 
consequence of increased local noise propagation in the model. 
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Figure A. 12: Mean environmental robustness Rj, against T for evolution of 
L = 10 targets in the Wagner threshold model for constant ICs. Both topologies 
demonstrate a significant increase from very low initial levels. The increase 
generated may be consequence of the selection for mutational robustness also 
affecting environmental robustness. 



Fundamentally the results are very similar to before with 
the Oikonomou and Cluzel ( 2006| l model. SFout topologies 
are consistently more evolvable than the ER topology, whilst 
both topologies increase their mutational robustness and envi- 
ronmental robustness under neutral mutation 



Appendix B. Calculation of order-chaos phase transition 

In order to calculate the order-chaos transition, we followed 
the method of Oikonomou and Cluzel] ( 2006| l. Let us consider 



two different sets of ICs, {(7^'} and {af^}- The difference be- 
tween these sets over time can be measured as with the Ham- 
ming distance as 



d(t) 



1 ^ 

N ^ 

1=1 



o'Pit)- 



ofm 



(B.l) 



The annealed approximation of IDerrida and Pomeau ( 19861 1 
measures whether noise propagation diminishes or extremities 
in the hmit of the Hamming distance between the two configu- 
rations tending to 



by selection. A significant increase in environmental robustness 
is observed in this case, although for both topologies, the typical 
maximum environmental robustness attainable starts levelling 
off around Rd ~ 0.6. This is different to the Oikonomou and 



|Cluzel| ( [2006| l model, where a higher environmental robustness 
was generated by both topologies. Again this difference may 
be a consequence of local noise propagation being stronger in 
the Wagner model. The inability to evolve more complete envi- 
ronmental robustness requires further research. The analysis of 
the attractor properties within the Wagner model would proba- 
bly be helpful. The result for constant ICs in Fig. A. 12 shows 



that under neutral mutation there is a clear optimisation for both 
topologies. The networks naturally possess very low environ- 
mental robustness atT - 1, particularly for the ER topology. 



M 



dd(t + 1) 



dd(t) 



limd(t)^0 



Y^Pikd.ki.pAki) (B.2) 



ki=\ 



where M measures the noise propagation in the network. When 
M > 1, the network is in the chaotic phase, whilst M < 1 is the 
ordered phase. For M - I, the critical phase is obtained. P(fe,) 
is the usual degree distribution, whilst Ps(ki) is the probability 
of local noise propagation. This measures the probability that if 
one of the inputs to a node is flipped at time t, this will change 
the output of the node at time f -i- 1 . For threshold models this 
is a function of kt, unlike in the Boolean function models. Here 
we calculate numerically the values of psik/) for both the update 
rules in the |Oikonomou and Cluzel| ( |2006| l model, as well as the 
Wagner model. These two probability functions are presented 
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Figure B.13: The probability of local noise propagation, Ps(k) against k for 
both the [Oikonomou and Cluzel (2006i model and the Wagner ( 1994 1 model. 
Local noise propagation is stronger in the Wagner model. Both distributions 
are approximately oc k^^-^ . 



Figure D.14: Mean mutational robustness R,,,, against T for evolution towards 
X = (10) with modified weight mutation method. As previously shown. SFout 
topologies possess greater robustness than their ER counterpart indicating that 
the previous mutation rales do not bias our result. 



Table B.6: The critical parameter values for the phase transition in the various 
topologies. 



Topology 



SF 

^^standard 



Oikonomou 



Wagner 



Kc = 3.54 
y, = 2.42 
Kr = 3.83 



Kc = 2.08 



in Fig. B.13 The Wagner model has much greater local noise 
propagation due to an "off" state still passing a signal to a node. 
Using these distributions the critical parameter values are cal- 
culated for the various topologies using Eq. ( |B.2| i and presented 
in Table Ij 



Appendix C. Total number of sequences for each repeat 
length 

There are a finite number of possible configurations a se- 
quence of bits of a given length may take. Two sequences are 
regarded as identical if they can be mapped on to each other 
through cyclic permutation. If a sequence can be broken down 
into a multiple number of shorter sequences, it is invalid. As 
such, the number of possible sequences A{n), for a sequence of 
a length n, is given by the following 



A{n) 



2"-T.keTkA{k) 



(C.l) 



Appendix D. Importance of mutation method with respect 
to increased robustness 

One possible explanation for the increase in mutational ro- 
bustness of the SFout topology would be a possible bias intro- 
duced by the weight mutation rules. This bias could arise in 



the following way. Given that any edge, incoming or outgoing, 
can have its weight mutated for the node under mutation, the 
probability of selecting an incoming node for a hub in the SFout 
topology will be much smaller Explicitly, there are many more 
outgoing edges that could be randomly selected for mutation. 
To assess whether this mutation method was introducing erro- 
neous results, we performed runs towards £, - {10) with modi- 
fied mutation rules. Instead of allowing the weight of any edge 
touching a given node to be mutated, we only allowed weight 
mutations of the edges incoming to that node. As both topolo- 
gies have the same incoming degree distribution, this removes 
any possible bias. 

A graph of the change in average mutational robustness, R,n, 
for 50 evolutionary runs is shown in Fig. |D.14| As can be 
seen there is limited difference between this and the original, 
demonstrating any bias does not affect the evolution of muta- 
tional robustness in this case. 



Appendix E. Population Diversity 

To determine whether the diversity of a population geno- 
types undergoing neutral mutations, we measured the genetic 
difference between every pair of inviduals each generation after 
maximal fitness is reached. The fractional genotypic distance 
D{A, B), between two networks A and B is given by 



D(A, B) 



1 
2M 



N N 

'=1 >=I 






(E.l) 



where T is the set of all factors of « and A(l) - 2 (i.e. and 1). where 



/w = 



[1 ifjc^O 
10 ifx = 



and M is equal to the average number of entries in a network's 
weight matrix, given by M = A^(^). Values for D(A,B) are 
taken over all pairs A, B e {population}, and then the values are 
binned in distances of 50/2M up to the maximal value of 1 . 
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Figure E. 1 5 : The heat of a square represents the fraction of pairs of networlcs in 
the population separated by a genetic distance of D at time T. As can be seen, 
the genetic distance gradually increases as neutral mutations are incurred to 
each network. The increasing "hot" lines demonstrate the population diffusing 
through genotype space creating a genetically diverse population. The sudden 
drops in diversity, occassionally observed, correspond to a new most recent 
common ancestor taking over the population. 



These measurements were performed on several different 
types of evolutionary run. A heat map of these distances is 
displayed in Fig. 



E.15 for a typical evolutionary run (once 



maximum fitness is reached) for evolution towards a X = {10} 
target. This figure demonstrates several features of an evolu- 
tionary run. The gradually increasing "hot" curves show the 
population spreading out from each other in genotype space as 
independent mutations are incurred. The gradient of the "hot" 
line decreases as mutations start occuring on top of each other. 
However, the "hot" lines can still increase to large values (tend- 
ing to D a; 1 for some lines), indicating that most of the muta- 
tions are neutral ones and lead to a large population diversity. 

An interesting feature of the plot are the occasional dramatic 
decreases in diversity. These are indicated through the collapse 
of a "hot" line to a much lower value on a new curve, corre- 
sponding to the introduction of a new most recent common an- 
cestor for the population. This ancestor may have been selected 
for or simply fixed through random sampling. 
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